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abstract 

The Galerkin finite element solutions for the 
scalar homogeneous Helmholtz equation are presented for 
no reflection, hard wall, and potential relief exit 
terminations with a variety of triangular element ori- 
entations. For this group of problems, the correlation 
between the accuracy of the solution and the orientation 
of the linear triangle is examined. Nonsymmetric ele- 
ment patterns are found to give generally poor results 
in the model problems investigated, particularly for 
cases where standing waves exist. For a fixed number 
of vertical elements, the results showed that symmetric 
element patterns give much better agreement with cor- 
responding exact analytical results. In laminated wave 
guide application, the symmetric pyramid pattern is 
convenient to use and is shown to give excellent 
results. 

NOMENCLATURE 

A constant coefficient 

Ag element area 

B constant coefficient 

c speed of propagation 

D domain of problem 

{F) global boundary condition vector 

H height of wave guide 

i -1 

K?®) local stiffness matrix component, Eq. (49) 

1 » J 

local stiffness matrix 
[K] global stiffness matrix 

k wave number, Eq. (2) 


L length of wave guide 

Nj local element shape factor 

[N^®^] local element shape matrix 
n total number of nodes 

n unit outward normal 

R residual error 

S boundary line of domain D 

t time 

global basis function Wi(x,y) 

X axial coordinate 

y transverse coordinate 

a/an derivative in normal direction 

7 ^ laplacian operator 

7 gradient operator 

Cg exit impedance, Eq. (11) 

♦ potential 

nodal value of potential 
value of potential in element, ()>^®^(x,y) 
( 4 .} global vector of nodal values 

element vector of nodal values 
u angular frequency 

Superscripts 

( ) approximate value of ( ) 


1 



(e) element value 


are made as to the advantages of the respective element 
orientations. 


dimensional quantity 


GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 


Subscripts 

e exit or area 

i index 


The governing differential equation describing 
the propagation of a harmonic scalar perturbation 
potential (|> is the classic Helmholtz equation. The 
two-dimensional form of the homogeneous Helmholtz equa- 
tion can be written as 


j i ndex 

INTRODUCTION 


3 ^ + 
ax 
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( 1 ) 


The Helmholtz equation models many physical 
systems. In acoustics, the classic wave equation 
reduces to the Helmholtz equation which describes the 
propagation of harmonic pressure disturbances (1). 
Electromagnetic propagation in wave guides is aTso mod- 
eled by the solution of the Helmholtz equation (2^). In 
fluid mechanics, a spectral approach to the Navier- 
Stokes equation also leads to the second order elliptic 
Helmholtz equation {3). The present paper is concerned 
with linear triangle finite element solutions to the 
Helmholtz equation. In particular, the optimum orien- 
tation of the triangle element pattern is examined for 
a number of typical boundary conditions associated with 
both acoustic and electromagnetic wave guides. 

In acoustics and electromagnetic finite element 
modeling of the Helmholtz equation, first order linear 
triangular elements made their initial appearance in the 
late 1960's and early 1970's. Recently, higher order 
elements have supplanted the first order elements 
because of their higher accuracy. Nevertheless, first 
order elements continue to find use because they are 
simpler to understand and program (less chance for ini- 
tial errors) and consequently they are an attractive 
starting point for a new research area. Many triangular 
patterns have appeared in the literature. In beam ana- 
lysis, for example, the nonsymmetric pattern of tri- 
angles shown in Fig. 1(a) was shown to give less error 
than the symmetric pattern shown in Fig. 1(b) (£ p. 71). 
At the present time, however, the literature does not 
suggest which triangular element pattern might best be 
employed in the solution of the Helmholtz equation with 
boundary conditions typical of duct propagation 
problems. 

Herein, the finite element solutions of the 
Helmholtz equation are obtained for plane wave propa- 
gation in a rectangular duct with a variety of linear 
triangle orientations. The correlation between the 
accuracy of the solution and the orientation of the 
linear triangles is examined for the following termi- 
nation conditions: (1) no reflected wave at exit; 

(2) hard wall (infinite impedance); (3) potential relief 
(zero potential). Boundary conditions (1) and (2) are 
Neumann conditions and enter the problem as natural 
boundary conditions. On the other hand, the potential 
relief constraint enters the problem as a Dirichlet or 
forced boundary condition. Conditions (2) and (3) give 
rise to resonance inside the duct. 

The first section of this paper presents the 
differential equation describing a scalar potential 
propagation in a wave guide. The second section pre- 
sents model problems along with closed form analytical 
solutions which will later be compared to the finite 
element results. In the third section, the finite ele- 
ment formulation of the model problem is given. In the 
fourth section, finite element solutions and analytical 
solutions are compared for a number of triangular ele- 
ment patterns and termination boundary conditions for 
the proposed model problems. Finally, recommendations 


where k is the wave number representing a ratio of 
frequency to a propagation speed. 

k=^ (2) 


where all the variables are assumed dimensionless. The 
nondimensional ization begins with the speed of propa- 
gation c*, frequency a and dimensional potential 
and introduces their nondimensional equivalents. 

The superscript ( ) denot|s a dimensional quantity. 

The speed of propagation c is normalized with 
respect to Cq the velocity of propagation in a 
referpncp medium. ]^he frequency u is normalized 
by H /Cq where H is a characteristic 
l|n^th. The potential is normalized with respect to 
CqH . In addition, lengths are scaled by H . 

In wave guide propagation problems, the boundary 
conditions will depend on the type of source forcing the 
propagating or standing wave, the termination of the 
guide, and the properties of the walls of the wave 
guide. Typical boundary conditions for acoustic propa- 
gation can be found in Ref. 5 while the boundary condi- 
tions in electromagnetic theory are discussed in Ref. 2 
(p.-67). The boundary conditions associated with the 
model problem under consideration in this paper will be 
presented in the next section. 

MODEL PROBLEM 


Finite element solutions of the Helmholtz equation 
will be determined for the two-dimensional boundary 
value problem shown in Fig. 2. In this problem, a uni- 
form potential is assumed at x = 0, 

♦(0) =1 (3) 


and the normal gradients of the potential are assumed 
zero at both the upper and lower- walls. 


3<t, 

1 ^ 1 

= 0 

an 

II 

O 

y=0 
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34> 

= 0 

an 

y= 

=1 ^ 

y=l 


(4) 

(5) 


At the exit (x = L) of the wave guide, the three 
separate termination conditions shown in Fig. 2 will be 
considered. Their significance will be discussed 
shortly when the closed form analytical results are 
developed. 


Closed Form Solution 

Since the boundary conditions are uniform at the 
entrance and exit of the wave guide and are identical 
along the top and bottom of the guide, the scalar (|> 


2 


will not depend on the y coordinate. Consequently, 
the governing Helmholtz equation reduces to 


along with condition (3) yield an equally simple 
solution 


d ♦ a. , 2 

— ^ + k A = 


dx^ 


k ♦ = 0 


which has a general solution of the form 




(6) 


(7) 


cos k(L - x) 
* ~ cos kL 


( 15 ) 


In this solution, both constants A and B have non- 
zero value. Consequently, the solution represented by 
Eq. (15) is. a combination of forward (e”^'^’') and 
backward (e^'^’^) traveling waves. A resonance rein- 
forcement of the waves occurs when 


If the time dependence of the. wave (neglected here) is 
assumed to be e^“^, the Ae“^'^^ represents a wave 
propagating in the positive x direction while 
BeiKx represent a wave moving in the negative x 
direction (exactly the opposite would occur had a 
e-iut dependance been assumed). Details on the 

nature of wave propagation can be found in any text on 
acoustics (6). 


kL = Y> • • • (16) 

The specific problem to be considered herein will 
again be for a k value of 2ir. In this case, the 
analytical solution becomes 

<t> = cos 2ir(L - x) (17) 


Exit Terminations 

Separate analytical and numerical solutions will 
be associated with each of the exit boundary conditions 
shown in Fig. 2. The significance of each of these 
conditions will now be discussed. 

(1) No reflections at exit (gp = 1) . Consider 
the case where waves only propagate to the right at 
X = 1. In this case, 

B = 0 (8) 

and the' exact analytical solution can be written as 
— i k X 

4 > = e = cos kx - i sin kx (9) 


In contrast to Eq. (13), ♦ does not have any imagi- 
nary component, so the real designation is not required. 

The term hard wall exit comes directly from acous- 
tics. For a hard exit wall, the acoustic velocity 
represented by the gradient of the potential will be 
zero. In this case, substituting Eq. (14) into the 
definition for the impedance given by Eq. (11) yields 

Cg = ” (18) 

Equation (18) will be conveniently employed in the 
finite element solutions to follow. 

(3) Relief Exit (» = 0) . The exit condition 

,(, = 0 at X = L (19) 


The coefficient A has obviously taken on a value of 
unity to satisfy the entrance condition (3). 

The exit termination boundary condition associated 
with this condition can be found by differentiating 
Eq. (9) and setting x = L: 



= -ike" 


■ ikL 


-ik4,(L) 


(10) 


In wave propagation problems, it is customary and 
convenient to define an impedance of the form 




or 




( 11 ) 


along with condition (3) yields a solution of the form 

♦ = k > 0 (20) 

Again, because both constants A and B in Eq. (7) 
have nonzero values, the solution given by Eq. (20) 
represents the sum of forward and backward traveling 
waves. In this case, the resonance condition occurs 
when 


kL = V, 2it, 3 ti, . . . (21) 

For this boundary condition, three specific example 
problems will be developed. For the example shown in 
Fig. 2, where L = 1 : 


For this definition, substituting Eq. (10) into 
Eq. (11) yields 


c 


e 


= 1 


(12) 


(a) k = 0 <, = 1 - X (22) 

(b) k = |- 4 i = sin |•(1 - x) (23) 


as shown as exit condition (1) in Fig. 2. 

The specific problem to be considered herein will 
be for a k- value of 2n. In this case, the real part 
of the analytical solution becomes 

Real ( ifi) = cos 2irx (13 ) 

(2) Hard wall exit a^i/ax = 0 (^ = °°) . 

The exit condition 



(14) 


(c)k=|i ♦ = sin _ x) (24) 

In a sense, boundary conditions (2) and (3) in 
Fig. 2 are quite similar in that they both lead to 
standing wave patterns which are identical over certain 
regions. However, boundary condition (3) was intro- 
duced because it enters the finite element solution as 
a forced (Dirichlet) condition in contrast to boundary 
conditions (1) and (2) which enter the problem as nat- 
ural boundary conditions. As pointed out clearly by 
Silvester and Ferrari (2, p. 11) these different formu- 
lations can significant^ effect the accuracy of the 
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finite element results. • In the' weak formulation (7_, p. 
140) of the finite element equations the order of the 
differential equation is reduced using integration by 
parts (Green's vector identity). In this case the nat- 
ural (Neumann) boundary conditions are not satisfied 
exactly but only in a certain mean-value sense which 
causes a contour integral in the finite element solution 
(to be presented shortly) to take on a prescribed value. 
The "strong" formulation, which enforces gradient 
boundary conditions exactly at a node ( 8 ), requires the 
use of Hermitian shape functions and is not considered 
herein. 

FINITE ELEMENT THEORY 

The finite element formulation is now generated by 
using the Galerkin method (9 to 11) to obtain an inte- 
gral form of Eq. (1) over tlTe whFTe (global) domain D 
shown in Fig. 2. A simplified and complete explanation 
of the Galerkin method is difficult to find in the lit- 
erature. The presentation here is given in fair detail 
and is tutorial in nature. 

System Discretization 

The continuous domain 0 is first divided into a 
number of discrete areas staked out by the nodal points 
as shown in Fig. 1(a). Although the nodal points are 
shown evenly spaced in Fig. 1(a), the advantage of 
finite element theory is that it allows the placement 
of the nodes at any position desired. Next, the con- 
tinuous potential 4 .(x,y) will be approximated (curve 
fitted) in terms of the nodal potential values 
located at x^,yj, see Fig. 1(b). The curve fit is 
required since we will integrate the governing Eq. (1) 
over the domain 0 shown in Fig. 2. This contrasts 
with the simplest form of finite difference theory 
(Taylor series) which usually only determines the 
potential at the lumped nodal points. 

Global Weighted Residual Approach 

In the classical weighted residual manner, the 
potential <)i is curve fitted by expanding in terms of 
the nodal values <tii(xi,yT) and a series of 
basis (shape) functions, such that 

n 

4>(x,y) = L Wi(x,y)*i = [W]{4,} (25) 

i=l ^ ’ 


where the basis or weight functions W-j(x,y) charac- 
terizes the spatial dependance of (f, and <(),• 
represents the unknown value of the potential ^ at 
specific nodal points in the global region. The hat 
over the .f indicates that it is the approximate 
solution to ♦. For example, iji-j would be asso- 
ciated with the potential value at the nodes shown in 
Fig. 1(a). In Fig. 1(a), the number of elements m has 
a value of 24, while the total number of nodes n has 
a value of 20. The nodes are numbered 1,2, . . . n 
and the global vector { 41 } represents the scalar values 
of the unknown potential 41 at each node, such that 
in matrix form 


(*1 = 


nj 

i>2 / 7 

or 4-2 


(26) 


The weight W-j has the property of being unity 
at node i and identical to zero at the other nodes. 

The weights are all assumed to be known func- 
tions; consequently, we must guess some approximate 
form of Wi before the weighted residual approach can 
be applied. The values of W^(x,y) at points other 
than the nodes will generally be finite nonzero values. 
However, in the finite element (Galerkin) approximation 
for W-j(x,y) to follow, throughout most of the domain 
the weight is assumed zero not only at all other 
nodes but also at all values of x and y outside the 
element under consideration. This Galerkin approxima- 
tion will considerably simplify the required 
integrations. 

In general, substitution of Eq. (25) into the gov- 
erning Eq. (1) and integrating over the domain D will 
not be equal to zero but will leave a residual error R. 




k^ 4 i] dx dy = R 


(27) 


In accordance with the method of weighted residuals, the 
assumed basis functions Wi and the distribution of 
errors R are forced to be orthogonal (R = 0) within 
the region, such that 

JJ" W^[ 7 ^ 4 i + k^ 4 i] dx dy = 0 

D 

yy* W 2 [v^ 4 > k^ 4 i] dx dy = 0 

0 


\ (28) 

jy W^[ 7 ^ 4 > + k%] dx dy = 0 I 

D I 


jy W^[ 7 ^ 4 > + k^*] dx dy = 0 I 

D 


Thus, there are n equations for the n nodal 414 
unknowns. In a direct analogy to the finite difference 
weighted residual control-volume formulation ( 12 , pg. 
30), each of the above equations can be thoughT”of as a 
higher order difference approximation at the nodal point 
where W, has a value of unity. Equations (28) can 
be written in compact form as 




dx dy = 0 


(i = 1, 2 . . ., n equations) (29) 


By making use of Green's vector identity (integra- 
tion by parts - 9, Eq. 9) and using the divergence the- 
orem of Gauss (^, p. 79, Eq. 4.7(b)), Eq. (29) becomes 

k^'^-j)dxdy - ^W-jv'T') *nds = 0 
D S 


(i = 1, 2 . . ., n equations) (30) 

where n is the unit outward normal and S is the line 
integral around the global domain D as shown in 
Fig. 2. 


4 



(33) 


In effect, we have reduced the order of the dif- 
ferential equation allowing us to employ the weak for- 
mulation of the finite element theory. 

Thus, in the Galerkin finite element approximation to 
follow, simple class Cq shape functions can now be 
used to approximate W-j . Across an element shown in 
Fig. 1 (a) for example, the class Cg functions are 
only continuous in the dependent variable ♦ and are 
discontinuous in slope. As previously mentioned in the 
discussion of the boundary conditions, if Eq. (29) is 
treated directly, then Hermitian functions are required 
to approximate Wi (8). In these functions both the 
variable and its slope are continuous across a boundary. 

Finite Element Approximation 


Both the specifications of the global weighting 
functions and the global integration over the 
whole domain D required by Eq. (30) are not practical. 
However, the integration can readily be performed by 
subdividing the domain into smaller elements Ae, 
defining the global shape function W-j in terms of 
the nodes of an individual element, integrating over an 
individual element and summing all the elements 
together. 

Equation (30) is valid over the entire domain D 
shown in Fig. 2 or any subdomain Ae, as represented 
by the area of a small triangular element embedded in 
the region as depicted in Fig. 2. To begin the finite 
element aspect of the weighted residual method, the 
domain D is assumed to be divided into m elements 
defined by n nodes, as shown by the example in 
Fig. 1(a). In this case, Eq. (30) can be rewritten as 



(i = 1, 2 . . ., n equations) (31) 


♦l 

*2 : 
*3 

and the shape function 


The form of the local 
on the type of element used. For the linear tp" angle 
element employed herein, the known value of 
is simple in form and can readily be found in nearly 
every text on finite elements (2^, Eq. 3.05). Like its 
global counterpart W-j, the local has the 

property of being unity at node j and zero at the 
other nodes in the triangle. For example, the magnitude 
of in element (1) is illustrated in Fig. 3(a). 

In this case, the magnitude of is represented 

by the height above the x-y plane. Np) for other 
selected values are also seen in Fig. 3. 

Notice that Nff' will have a nonzero value for 
elements (12), (13), (14), (19). (20), and (21). A 
superscript is required for to denote which 

element it resides, since the magnitude of this values 
the shape function at any x and y inside the ele- 
ment will depend on the location of the other corner 
nodes. 

Replacing 'J' by in Eq. (31), the new 

governing equation becomes 



{*} 


( 20 ) 


1*13 

^7 

hs' 




- LN 13 N 17 N 18 J 


(34) 


shape matrix [n(®^] depends 



Local Shape Factors 


[i = 1, 2, 


n equations) 


(35) 


Each element is defined by the nodes around its 
perimeter. To represent the variation inside the ele- 
ments of the field variable ♦ and its derivatives, 
local interpolation shape functions N-j(x,y) for the 
linear triangle are written in both scalar and matrix 
form as 


(|>(x,y) = 4 i^®^(x,y) = N[®^(x,y)(|i|®^ + N^®^(x,y)*^®^ 


- N(^)(x,y)4^) 

j=l 




(32) 


where is the vector of nodal values of 

^ for a general element e with subscripts 1, 2, and 
3 representing the nodal positions, as shown in Fig. 2. 

For example, for element 20 in Fig. 1(a), the vec- 
tor would take on the form 


Galerkin Approximation 


In general, the approximation for the dis- 
tribution and the weight W-i can be different, as 
seen by now comparing Eqs. (25), (32), and (35) and, for 
example, as usually applied in weighted residual con- 
trol-volume finite difference theory (J^, pg. 32). In 
conventional finite element analysis as applied here, 
however, the weight function and the profile function 
are usually equated. The weight W-j(x,y) is now 
approximated by multiple values of Adapting 

the Galerkin approximation to the more general weighted 
residual approach assumes that 

Wj(x,y) = Nj®^x,y) (36) 

in all elements containing the node i. Approximation 
for Wi can be visualized by the pyramid-like shapes 
displayed in Fig. 3. For all elements which do not 
contain the node i, the weight Wi is assumed zero 
not only at all other nodes (as required by the general 
definition of Wi ) but also at all values of x and 

y- , , 

Recognizing that N)®) is zero for all ele- 
ments not having the unknown <>.j associated with a 
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particular element, as shown in Fig. 3, the finite ele- 
ment Eq. (35) can be written in compact form as 



( i = 1, 2, . . n equations) (37) 

where N-j are the known shape factors (^, Eq. (3.05)) 
Surface Integral 


The surface integral in Eq. (30) can be rewritten 
along the boundary in Fig. 2 as 

iw.V4,*nds = - y* W. ||dy - /w. 

S FRONT ' BOT ^ 


* y*w. i^y + y*w. 

exit TOP ^ 


3y 


(38) 


Since, 34>/3y is zero (Eqs. (4) and (5)) along the 

top and bottom for our 

model problem, Eq. (38) simplifies to 


X.v**nds = - X W. dy + /w,. 

S ’ FRONT ^ exit ’ 


dx 


dy 


(39) 


or, utilizing the definition of impedance given by Eq. 
(11) yields 


L 


..V4i*nds = - 
S FRONT 


f ^ _lk r 


W.^dy 


(40) 


exit 


Employing 4i(®) and the Galerkin approximation, 
Eq. (36), Eq. (40) becomes 



- (41) 

® exi t 

In terms of the general in Eq. (40), both 
surface integrals would contribute to every element in 
the domain. However, with the Galerkin approximation 
the first term on the right hand side of Eq. (41) will 
only contribute when the nodal point i is on the right 
hand face. In Fig. 1(a), this would be points 1,2, 3, 4, 
and 5. For all other nodal positions, the value of 
Ni on the front face will be identical to zero by 
nature of the Galevkin approximation to shown in 
Fig. 3. However, at points' i on the front face 
boundary, (|> (Dirichlet) boundary conditions are 
specified (forced). Thus, the i'"" equation of Eq. 

(37) representing a boundary point will be discarded. 

The exact known answer will take its place; 

-■ = 1 (42) 


Thus, without loss of generality the first term or the 
right in Eq. (41) can be discarded, and 

= / nS"^("V (43) 

S exit 


For the relief exit condition, Eq. (19), the surface 
integral at the exit would also be discarded for the 
same reason. 

The finite element equation for a general i node 
can now be written as 


m 






® exit 


(i = 1, 2, . . ., n equations) (44) 


EXAMPLE 

The original global weighting function W-j 
associated with the i^^ global node has now been 

assumed to be composed of the separate N\ 'shape 
functions which reside in the local elements surrounding 
the i^^ node, as shown in Fig. 3. For example, for 
= N^^^, only the (5), (6), and (7) elements will 

contribute to the evaluation of the summation for the 
fourth node. Thus, for the grid and element notation 
of Fig. 1(a), the general Eq. (44) can be expanded into 
the form; 


- 1 

i=l 

4.2 = 1 

i=2 

*3 = ^ 

i=3 

4.4 = 1 

i=4 

♦5 = ^ 

i=5 


L 

e=l,2,9 



- k^*^®^N^®^)dA = 0 i=6 


e=2,3,4, JJ'^ 
9,10,11 Ae 




>2:(e)M(e) 


)dA = 0 i=7 


(equation continued on next page) 


6 


24 


e=22 


fi 

Ae 




. JJS = 0 'i=19 


II 

Ae 


Se 




*r/' 


n ( 24 );( 24 ),^ ^ 0 .^20 


Se 


(45) 


In this example, the first five equations asso- 
ciated with Eq. (44) reduce to the forced Dirichlet 
boundary condition. The next two equations are typical 
of Eq. (44) applied to the central regions. No surface 
integrals appear since the, value of. N, representing 
in the domain D is zero outside the element 
containing i along the exit surface. Finally, when 
the i node coincides with the exit boundary, a sur- 
face integral .contribution occurs, as in the last two 
equations. For i = 20, no summation is required since 
only one element area is associted with the final node. 
There are 20 equations associated with the 20 unknown 
values of the potential. Obviously because of the known 
boundary conditions, the equations could be reduced to 
fifteen. 


Finite Element Equations 

Expanding in terms of the weighting func- 
tions, Eq. (32) and neglecting the surface integrals for 
simplicity, -Eq. (44) becomes 


^ /■ r 

- k^N^^^N^^hjdAg = 0 

(i =1, 2,' . . ., n equations) (46) 

where the local potential vector can be 

pulled outside the area integration since it does not 
depend on the area coordinates. 

. Equation (46) could in principle be evaluated for 
each i to yield a set of n simultaneous equations 
for the n unknown values of (f-j, as in Eq. (45). 
However, this approach is not readily implemented on a 
digital computer. Rather, each element is treated 
independently and the contributions of a single element 
to all the equations in Eq. (46) are determined simul- 
taneously from the following: 



[vN^^h 




dA = 0 


(j = 1, 2, 3) 


(47) 


In contrast to the i subscript of Eq. (46) for all the 
nodes in the global domain, here the subscript j sums 
only the three nodes of a particular element. In com- 
pact standard finite element form, Eq. (47) can be 
written as 




= 0 


where 


4!1=/A ["S'M:' 






dx dy 


(48) 


(49) 


For triangular elements, the evaluation of Eq. (49) is 
quite simple and presented in many texts (11, p. 149, 
Eq. 8.42 for the first two terms and p. 45Tbr the last 
term). The surface integral is evaluated in a similar 
manner. 

Using the standard finite element procedures (11), 
the elemental Eqs. (48) can be assembled into Eq. (TTB") 
and the boundary conditions applied to yield the 
following global set of simultaneous equations which we 
can solve to obtain the unknown nodal potentials: 


[K]l4i} = IF) 


(50) 


Here, the column vector IF) contains boundary con- 
dition information and the global stiffness matrix [K] 
is the sum of the known local stiffness matrices 
[k(®)]. The solution of Eq. (50) by a banded Gauss 
solver yields (j>j. Except for a few minor modifi- 
cations, a complete fortran listing of the banded solver 
can be found in Ref. 13, program 16, page 64. 

Next, the solutions of Eq. (50) for a number of 
element discretization patterns and boundary conditions 
will be compared to the analytical solutions. 

RESULTS AND COMPARISONS 


Finite element solutions for a variety of exit 
terminations will now be compared to the closed form 
analytical solutions presented earlier. Four distinct 
triangular grid orientations will be considered. 

Nonsymmetric Grid 


The nonsymmetric grid shown in Fig. 4 was used to 
calculate wave propagation without reflection at the 
exit, Eq. (13). This grid orientation is labeled hon- 
symmetric since views looking from the top and bottom 
will be different. In this example, a rule of thumb [5) 
requires that approximately 12 linear triangles be 
employed to accurately resolve one complete harmonic 
oscillation. In the vertical direction, only four tri- 
angles were chosen since no variation of ((> in the 
y direction was expected. 

Figure 5 displays a comparison between the finite 
element calculations, shown, by the solid symbols and the 
analytical calculation shown by the solid line. In this 
and all comparison figures to follow, the finite element 
calculation are displayed for the nodal points at upper 
and lower walls and at the centerline. The heading on 
the figures denote the respective positions. 

Reasonable agreement is seen between the finite 
element results and the analytical curves. However, 
errors as high as + 10 percent occur at selected posi- 
tions. Other gricTorientation to be presented will 
improve on this comparison. The deviation between the 
upper and lower nodes is a direct result of the non- 
symmetric triangle pattern. A different arrangement of 
nodes contribute to the global "difference" equation 
used to calculate at the upper and lower sur- 
faces for the same axial position.. 

Figure 6 displays similar results for the non- 
symmetric grid for the hard wall exit conditions 
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(3(j)/3x = 0). In this case, the errors are very large 
and the result not satisfactory. 

The results for the nonsymmetric grid for the 
potential relief exit condition (♦ = 0) are now con- 
sidered. For a wave number of zero (k = 0), exact 
agreement is obtained. The upper, center, and lower 
surface finite element results coincide because of the 
linear nature of the solution, Eq. (22). For k 
greater than zero (k = u/2), the results are very good 
and just a slight variation between the upper and lower 
surfaces are seen. For k = 3n/2 shown in Fig. 7, the 
results are totally unsatisfactory. 

In order to check for convergence of the non- 
symmetric grid when k = 3ir/2, the vertical direction 
was divided into 20 elements. For this case as shown 
in Fig. 8, the agreement between the finite element and 
analytical results is good. However, the use of many 
elements to resolve a vertically nonvarying potential 
is quite unacceptable. Therefore, let us now consider 
some symmetric triangular grids. 

Symmetric Diamond Pattern 

The symmetric diamond pattern shown in Fig. 9 was 
used to calculate the potential relief boundary condi- 
tion (41 = 0, k = 3n/2) as shown in Fig. 10 and the 
wave propagation without reflection at the exit as dis- 
played in Fig. 11. In Fig. 10, agreement between the 
finite element calculations and the analytical results 
are excellent, and thereby preferred over the non- 
symmetric grid displayed in Fig. 4. Only the upper wall 
values are seen, since they are identical to the lower 
wall values and thus obscure the lower values. In Fig. 
11 for the nonreflecting exit, some significant error 
is still seen at an x value of 0.5. 

In practical applications, wave guides might have 
both vertical and horizontal laminated changes in mate- 
rial. To more conveniently resolve the boundary between 
laminated regions, some additional symmetric triangular 
patterns will now be considered. 

Symmetric Saw Tooth Pattern 

The symmetric saw tooth pattern shown in Fig. 12 
was used to calculate the potential relief boundary 
condition (4. = 0, k = 3u/2) as shown in Fig. 13. As 
seen in Fig. 13, excellent results are obtained. There 
is however, a slight variation between the center values 
and the wall values. Again, only the upper wall values 
are seen, since they are identical to the lower wall 
values and thus obscure the lower value. 

Symmetric Pyramid Pattern 

The symmetric pyramid pattern shown in Fig. 14 was 
used to calculate the potential relief boundary condi- 
tion (4, = 0, k = 3ii/2) as shown in Fig. 15 and the 
wave propagation without reflection at the exit as dis- 
played in Fig. 16. In both cases, agreement between the 
finite element calculations and the analytical results 
are excellent. In these cases, the values of the 
potential for the upper and lower walls and centerline 
are identical.. Thus, only the upper wall data point 
appears in these figures. The agreement in Fig. 16 at 
X = 0.5 represents a significant improvement over the 
diamond pattern results of Fig. 11. 


CONCLUDING REMARKS 

The finite element solutions for the Helmholtz 
equation were presented for no reflection, hard wall, 
and potential relief exit terminations with a variety 
of triangular element orientations. Nonsymmetric ele- 
ment patterns were found to give generally poor results 
in the model problems investigated. For a fixed number 
of vertical elements, the results showed that symmetric 
element patterns give much better agreement with cor- 
responding exact analytical results. For layered grids 
which may be needed for laminated wave guide applica- 
tions, the symmetric pyramid pattern was found to give 
the best results and was very convenient to program. 
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(b) Symmetric x-y directions. 


Figure 1. - Linear triangular elements. 
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Figure 2 . - Mixed boundary value problem. 
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(b) N2 Shape functions \^N2 ,N2 .N2 j 

Figure 3 . - Finite element approximations of 
Global basis functions. 
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Figure 5.- Potential 0 profile using non-symmetric 
discretization for k= 27T and AX = 0,083 with non- 
reflecting exit, Ze = 1 and 5 vertical nodes. 
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Figure 6. -Potential 0 profile using non-symmetric 
discretization for k = Ztt and AX = 0. 083 with 5 verti- 
cal nodes and hard wall exit (d0/dx = 0. 




Figure 7.- Potential cp profile using non- 
symmetric discretization for dimensionless 
wave number k= 37r/2with 5 vertical nodes 
and Ax = 0. 1 with 0 = 0 at exit. 



1.6 


1.2 


NODAL FINITE ELE/VIENT VALUES 
I UPPER BOUNDARY y= 1.0 



THEORY, 
EO. (24) 
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Figure 8.- Potential (D profile using non- 
symmetric discretization for dimensionless 
wave number k ■ 37rf2 with 20 vertical nodes, 
Ax = 0. 1 and 0 “ 0.0 at exit. 



Figure 9. - Discretization of solution 
domain using symmetric diamond 
pattern for k = 27rand Ax = 0. 083 
with 5 vertical nodes (symmetric 
x-y direction). 



Figure 10. - Potential 0 profile using sym- 
metric diamond discretization for dimension- 
less wave number k = 37r/2with 4 verti- 
cal nodes, Ax = 0. 1 and 0 = 0 at exit. 
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Figure 11. - Potential 0 profile using symmetric diamond 
discretization for k = Zjrand Ax = 0.083 with non- 
reflecting exit, = 1. 




Figure 12. - Discretization of solution 
domain using symmetric saw tooth 
pattern for k = 3 ti/ 2 and Ax = 0. 1 
with 5 vertical nodes. (Symmetric 
y-direction. ) 



Figure 13. - Potential 0 profile using saw tooth 
discretization for dimensionless wave number 
k = 3 tt/ 2 with 5 vertical nodes, Ax = 0. 1 and 
(p = Oat exit. 




Figure 14. - Discretization of solution 
domain using symmetric pyramid 
pattern for k = 3ir/2 and Ax = 0. 1 
with 5 vertical nodes. (Symmetric 
y-direction. ). 
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Figure 15. - Potential (t) profile using symmet- 
tric pyramid discretization for dimensionless 
wave number k = 3 ti/ 2 with 5 vertical nodes, 
Ax = 0.1 and (P = 0 at exit. 
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Figure 16, - Potential 0 profile using symmetric pyramid 
discretization for k =2ir and Ax =0.083 with non- 
reflecting exit, = 1 and 5 vertical nodes. 
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